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Abstract 

We suggest a unified spectrally matched optimal grid approach for finite-difference and 
finite-element approximation of the PML. The new approach allows to combine optimal 
discrete absorption for both evanescent and propagative waves. 



1 Introduction 

We approximate the Neumann-to-Dirichlet (NtD) map of wave problem in unbounded 
domain. After Fourier transform we obtain, 

u xx — \u = 0, x G [0, oo] (1) 

and due to the infinity condition we are limited to outgoing wave solutions, 

u = ce~^ x 

satisfying NtD condition, 

u , 1 



u. 



l^o - (2) 

Here A = k 2 — u 2 , where k and u are respectively (tangential) spacial and temporal 
frequencies. Also, ([I]) can be equivalently rewritten in the first order form as, 

u x = sv, v x = su, (3) 

where s = y/\. In terms of u and v, condition (|2j) can be equivalently rewritten as, 

-U=o = -l. (4) 

v 

The NtD can be numerically realized via rational approximation theory using 
several approaches HOI HH [HI HI H21 HS1 ?, etc]. In [H [I] and pEJ [EE] this 
approximant was realized as respectively finite-difference (FD) and finite-element 
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(FE) discretization of an absorbing layer similar to well known Perfectly Matched 
Layer (PML) [BJ. In particular, the FD scheme was designed as an optimal rational 
approximant separately for evanescent solutions corresponding to A > [Hj and 
propagative waves [1] corresponding to A < 0, but not for the both types of the 
solutions simultaneously. On the other hand, FE approach is more flexible; while 
[T2] focuses on propagative waves, it was shown in [18] that both propagative and 
evanescent waves can be treated simultaneously. Most recently, these FD and FE 
approximations are interpreted as special quadrature rules with complete wavefield 
approximation... [?]. 

In this paper, we show that simultaneous treatment of propagative and evanescent 
waves is possible not only in FE setting, but also in FD setting. The key to this obser- 
vation is a recently-discovered equivalence between the FE and FD approaches for the 
two-sided problem. Utilizing this link, we present two alternative approaches to im- 
plement the NtD map and comment on their relative merits. Furthermore, utilizing 
Zolotorev approximation theory and complete wavefield approximation interpreta- 
tion, we present an NtD map that is an optimal approximation for propagating as 
well as evanescent waves. 

The outline of the paper is as follows. We start in section 2 with the overview of 
optimal rational approximation of the NtD map by considering both propagative and 
evanescent waves. Section 3 contains the description of FE and FD approximations 
of two-sided problems and the equivalence between them. In section 4, we consider 
rational approximation of the NtD map of the exterior problem and present FE and 
FD realizations. The implementation details in time domain and relative merits 
of the (or three) approaches are considered in section 5. Numerical examples are 
presented in section 6. Finally, section 7 concludes the paper with some closing 
remarks, (where do we fit complete wavefield approximation of Tom?) 



2 Optimal Rational Approximation of NtD map for Propagative 
and Evanescent Waves 

Let us for simplicity consider time-harmonic case with u — 1 and consider time- 
dependent problems later. Let us present our rational approximant of —A" 1 / 2 as 

-A- 1 / 2 ^ J R(A)=p(A)/g(A), (5) 

where p and q are polynomials of degrees K — 1 and K respectively. Introducing a 
polynomial of degree N = 2K given by, 

h(s) = sp(s 2 ) + g(s 2 ), 

with s = we transform ([HJ) to Newman function 

W = pM/«M = = I'tX't («) 

s[h(s) + h{— s)\ s[l + h(s)/h{— s)\ 

Then the relative error of the NtD map is approximately proportional to the reflection 
coefficient, 

h(s) 

h(- s y 
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According to (jlj), the exact solution 
due to the approximation error, 



In reality, 

h(s) 



s) 



i.e. 



(u, v) of (]3]) is proportional to (1, — 1} 
(u,u)|o = ci(l,-l) + c 2 (l,l), and a =_ 
the reflection coefficient is the ratio of the incoming and outgoing waves. 

Minimization of on a real positive interval is the classical first Zolotarev 



h{-s) 



problem solved in 1872. Zolotarev's solution was first applied to the optimal FD ap- 
proximation of the NtD map for evanescent solutions in [14] and then to the approx- 
imation of propagative modes in [1] . The ABC for both propagative and evanescent 
waves should approximate the true NtD map on both negative [—1, Ai] and positive 
[A 2 ,A 3 ] intervals. They respectively correspond to intervals S p = \/Xi\ and 

[S e = \f\~2, a/A~3] of variable s. The so called spectrally matched finite-difference 
scheme (a.k.a FD Gaussian spectral rule or optimal FD grid) [7J etc] allows arbi- 
trary h(s), but does not simultaneously treat propagative and evanescent waves. On 
the other hand, propagative and evanescent waves have been simultaneously treated 
using FE approximation in [18]. The specific approximation is based on linear FE 
approximation with midpoint integration [12], which is linked to special rational ap- 
proximation [II] with 

h(a) = t(sf, (7) 

where t is a polynomial of degree fcQ Hence, considering the success in JT8] , we 
limit the current treatment to the restricted form of h in ([7J). With such restriction, 

h(s) 



minimization of max 



H-s) 



is equivalent to solving 



mm max 



t(s) 



t-s) 



(8) 



It is well known that the necessary and sufficient conditions for optimality of a 
real rational approximant on a real interval is so-called the Equal Ripple Theorem 
(ERT) [in]. It says that the optimal error of [(K-l)/K] approximant has 2k — 1 
zeros and 2k equal absolute value alternating extrema on the interval of optimality. 
Generally, there is no similar result for complex rational approximation [17]. Here, 
instead of minimizing (jSJ), we construct an approximant based on classical Zolotarev 
results. We hope that its error is close to®. 

If t = t e t p , degt e = I < k, degt p = k — I, where t e and t p have respectively (non- 
coinciding) roots on S e and S p , then has 2k + 1 maxima on S p U S e . Moreover, 



te{s) 



te(-s) 



1 on Sp and 



t p (-s) 



t(-s) 

1 on S e , which implies that 



and 



max 



max 



t(s) 


= max 

Se 


te(s) 


t{-8) 


te(s) 


t(s) 


= max 

Sp 


t p (s) 


t(-s) 


t p (-s) 



1 However, as it will be shown in the Section 4, simultaneous treatment of propagating and 
evanescent waves is even possible with more general h(s), if the FD approach is used. 
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Thus, we can take as t e and t p as the classical optimal Zolotarev approximants on 
S e and S p respectively, and obtain the quality of the total approximation the same 
as the one of the separate problems. 

The remaining question is: can the constructed approximant be optimal in global 
sense, or, at least, how close is its error to (JS}. Obviously, max Sp 



t(s) 



*(«) 



*(-*) 



*(-«) 
and max^ 



and max5, 

(5) 



t( S ) 



t(-s) 



t(-s) 

for a count- 



may be different. Varying / one can equate maxs p 

able set of arrays Ai, A 2 , A 3 . 

Here we conjecture that the ERT can be extended to the first Zolotarev problem 
on two intervals in C in the following way. 

Conjecture 1: Let t e (s)/t e (—s) and t p (s)/t p (—s) be the solutions of the Zolotarev 
problems on S e and S p respectively. 

1. There are infinitely many arrays Ai, A 2 , A3 for which there exists I, such that 



max 

Se 



Us) 



max 



tp{s) 



(9) 



2. If © is valid, then t 



t e t p solves 



Results of [?] indicate that, if (Q is valid, then at least the approximant is optimal 
in the Cauchy-Hadamard sense. Generally, it is always possible to find / such that 



t(s) 



t(-s) 



and maxs e 



t(s) 



t(-s) 



are of the same order, in which case, it is natural to 



assume that the approximation error will be of the order of (jBJ). 



3 Equivalence of FE and FD Approximations for Two-sided 
Problems 

While the emphasis of this paper is on the approximation of the one-sided problem 
on [0, 00), in this section, we consider the two-sided problem on [0, 1] and show that 
there exist equivalence between spectrally matched FD grids and midpoint integrated 
linear FE mesh. We then utilize these results in Section 4 to construct an effective 
NtD map for the one-sided problem on [0, 00). 



3.1 Continuum problem 

QUESTION: You have used * for many row vectors and matrices. Should 
we be just using transpose? decide later. 

Let us consider eq. (jBD on [0, 1], and define the two-sided DtN map as matrix- 
valued function F(s) G C 2x2 

F(s)u b = v b , 

where u b = [«(0), u(l)]*, v b = [v(0),—v(l)]*. It is easy to see that (u,v) is a linear 
combination of, 

(e ±sx ,±e ±sx ), 
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and simple computation shows that 
1 



F s) 



cosh(s) — 1 
— 1 cosh(s) 



sinh(Y 

where Z is an orthogonal matrix 



tanh(s/2) 

coth(s/2) 



(10) 



V2 



1 -1 
1 1 



Similarly, we define propagator operator from left to right as matrix-valued function 
G(s) G C 2x2 G(s)w(0) = w(l), where w = (u,v)* and from (JTUJ) we obtain 



G 



cosh(s) sinh(s) 
sinh(s) cosh(s) 



exp(s) 
exp(— s) 



Z*. 



3.2 Discrete problem: linear FE mesh with midpoint rule 



It was shown in [12] that the discretization of the original second-order from in ([I]) 
with midpoint-integrated linear FE mesh would lead to exponential convergence of 
the NtD map. Furthermore, it was shown in [T3] that such a FE discretization is 
equivalent to Crank-Nicholson discretization of the first order form ([3]), i.e. 



U i+ i -Uj _ V i+1 + Vi V i+ i - Vi _ U i+ i + 

I ~ S 2 ' L " S 2 ' 



i = 1, . . . n. 



:i2i 



where / i; i = 1, . . . , n are the FE lengths with YH=\ h = 1- ^ can be easily verified 
that (uj, Vj), j = 1, . . . , n, is a linear combination of 

yr 1 ± ks/2 Jj l±l iS /2 
l\lTks/2' l\lThs/2 / 

Comparing the above (approximate) solution with the exact solution and noting that 
Y^i=i h = 1) ^he FE solution approximates the exponential as 



where 



Assuming 



exp(s) ~ exp(s) = t(—s)/t(s), 

n 

^(^)=n( 1 -^/ 2 )' 



i=l 



u(0) = Ul, u(l) = u n , v(0) W Vi, v(l) w v n 
we can compute the approximate NtD map as, 

1 



F(s) 



sinh(s) 



cosh(s) —1 
— 1 cosh(s) 



tanh(s/2) 
coth(s/2) 



(13) 



Z*. (14) 
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Here, 
sinh(s) 



exp(s) — exp(-s) 



tanh(s/2) 



t(s)-t(-s) 



• 1/ \ , / x exp(s) + exp(-s) 
sinn(s), cosli{s) = ~ cosh(sj 



tanh(s/2), coth(s) = l/tanh(s) ~ coth(s). 



t(s)+t(-s) 

Similarly, the discrete propagator from left to right matrix can be computed as 



G 



cosh(s) sinh(s) 
sinh(s) cosh(s) 



exp(s) 
exp(-s) 



(15) 



Vectors -J=(l,±l) are the eigenvectors of G, so it has so called fixed point property, 
i.e., if u(0)/v(0) = ±1 then u(l)/v(l) = ±1 and vice versa. This implies that, if 
exact half-space BC @ is applied at x — 0, it will be also valid at x — 1 regardless 
of the accuracy of the FE approximation. In other words, adding an FE-discretized 
interval to a half-space does not alter the NtD map of the half-space. Furthermore, it 
was shown in [12] that adding a midpoint-integrated finite element to an approximate 
half-space can only decrease the approximation error in the NtD map. This property 
was used in [18] to enhance the approximation, originally designed for propagative 
waves, to simultaneously absorb evanescent waves. 



3.3 Discrete Problem: spectrally matched finite-difference grids 

It was shown in [] that one-sided, two-point BVP can be solved with staggered FD 
method with exponential convergence at the end points. The main idea was to link 
the staggered FD approximation to rational approximation of the exact NtD map 
and optimizing the resulting approximation using Zolotorev theory. This method 
was later extended to the solution of the two-sided problems by splitting the solution 
into odd and even parts and solving two one-sided problems on half-intervals using 
dual grids. Formerly called optimal FD grids, the basic idea of spectrally matched 
FD grids is summarized below. 

Let us introduce the FD grid steps hi, hi, i — 1, . . . , k. We split the DtN map into 
odd and even parts and compute each of them using a FD scheme on half interval. 
The odd and even problems can respectively be written in mutually dual form as: 
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It is known [7] that, 



Ui 



h\s 



his + 



h 2 s + ... 



hk-is + 



h k s + 



h k s 



Combining odd and even parts we obtain, 
u(0) = ul + u°, u(l) 



Ui 



u 



v(0)+v(l) 



v{0) -v{l) 



i- 



and the FD-NtD as 



fk 







l/fk 



z*. 



(17) 



(18) 



(19) 



Construction of spectrally matched grids involves a reverse procedure. First, 
rational approximation theory is used to obtain fk that approximates the NtD map. 
The resulting rational function is then used in (I17j) to compute the grid steps hi, hi 
using simple **** algorithm []. **** algorithm also constructively shows that any 
[2k-l/2k] rational function can be converted into an equivalent FD grid. 



3.4 Equivalence of discrete problems 

In this section, we show that if the number of finite elements are chosen to be even 
(n = 2k), the approximate NtD maps from FE and FD grids are equivalent. 

Lemma 1: For any set of parameters U G C, I = l,...,2k there exist parameters 
hi, hi G C U oo, I = 1, . . . , k, such that 

f(s) = tanh{s/2) (20) 

and vice versa. 

Proof. For any set of parameters hi, hi G C U oo there exist polynomials p and q (at 
most) degree k — 1 and k respectively, such that the continued fraction expansion 
ffTTj) can be presented as 

_ sp(s 2 ) 
h q(s 2 ) ' 

and vice versa. Equating numerator and denominator of fk and tanh we equivalently 
transform ( 120]) to polynomial identities 

sp(s 2 ) = t(s) - t(-s), q(s 2 ) = t(s) + t(-s). (21) 

Since p, q and t can be arbitrary polynomials of degree k — 1, k and 2k respectively, 
then for any p, q there is t satisfying (1211) and vice versa. (Aren't these different p 
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and q? If so, it is important not to confuse the polynomials p and q with 
the polynomials in the second equation in section 2. Should be rename 
these as p and g?) Agree, will do it later □ 
From the lemma, ( TH|) and ffT9l) . we obtain the following result about equivalence 
of the FE and FD DtN maps. 

Proposition 1: If ( 120]) is valid, then 

F(s) = F(s). 

Formula (|2"T|) can be used for computing the equivalent FE from the FD and vice 
versa. 

If the DtN maps are identical, then formula ()15p can also be used for computing 
the propagator matrix for the FD approximation. 

If exp(s) matches exp(s) in n non-coinciding frequencies, then fk matches tanh(s) 
at the same frequencies, and fk is Stieltjes function, hi, hi are real positive, and the 
problem becomes Hermitian. 



4 Approximation of exterior problems 

Let a discretized interval Qi = [a?_,x + ] have the propagator matrix (from left to 
right), 

expi(s) 
expi(-s) 

in the spectral coordinates, where expi(s) = ti(—s)/ti(s) defined as in the pre- 
vious section. First, let us impose the reflection coefficient h 2 (s)/h 2 (—s) at x + , 
i.e., at the right boundary any nontrivial solution can be represented as w(x + ) = 
[ch 2 (s), — ch 2 (— s)]* in the spectral coordinates, where c ^ is an arbitrary constant. 
Then the reflection coefficient at the left boundary will be the ratio of the components 
of = Q~ 1 w(x + ). That is, 

V ' h 2 {s) hi-s) 2 h 2 (-s) v ; 

If we impose the Dirichlet condition at the right boundary of Qi, which corresponds 
to h 2 (s) = 1, then the reflection coefficient will be f ^^ 2 . Let us now assume that 
we have a connected interval Q = Q\ U Q 2 with the Dirichlet condition at the right 
boundary (£l 2 is assumed to be on the right), and /^(-l) ^ s ^ ne reflection coefficient 
of n 2 . Then ( |22l) would yield the reflection coefficient of Q that is just the product 
of the reflection coefficients of the two subdomains. 

Now, let as assume, that we use the discrete problem in Q for the approximation 
of (jlj), i.e., h(s) from fl6]) can be presented as h(s) = ti(s) 2 h 2 (s). If we set h\ = t\ 
and t± = t e and t 2 = t p , then the reflection coefficient of Q will be identical to the 
one discussed in Section 2. However, Dirichlet condition on the right of Q 2 makes it a 
one-sided problem, and it is not necessary to restrict to the two-sided approximation 
in the previous section. In fact, the original FD optimal grids are optimized for 
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the one-sided problems and can be used effectively for f2 2 . This approximation is 
equivalent to the odd part of the FD approximation (I16p . i.e., 

= stt,, 2 = 1, . . . , k, u k+ i = 0. 



io onnalihr r, I e 1 — _ 



Then /i 2 can be obtained from the equality /&(s) = r4^rr4z4 ' ^ e -' ^ can ^ e an 



arbitrary polynomial of degree 2k. { 
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